home *** CD-ROM | disk | FTP | other *** search
/ Chip 2007 January, February, March & April / Chip-Cover-CD-2007-02.iso / Pakiet bezpieczenstwa / mini Pentoo LiveCD 2006.1 / mpentoo-2006.1.iso / livecd.squashfs / usr / lib / python2.4 / site-packages / Numeric / LinearAlgebra.py < prev    next >
Text File  |  2006-01-20  |  18KB  |  516 lines
  1. # This module is a lite version of LinAlg.py module which contains
  2. # high-level Python interface to the LAPACK library.  The lite versioho
  3. # only accesses the following LAPACK functions: dgesv, zgesv, dgeev,
  4. # zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf, dpotrf.
  5.  
  6. import Numeric
  7. import copy
  8. import lapack_lite
  9. import math
  10. import MLab
  11. import multiarray
  12.  
  13. # Error object
  14. class LinAlgError(Exception):
  15.     pass
  16.  
  17. # Helper routines
  18. _lapack_type = {'f': 0, 'd': 1, 'F': 2, 'D': 3}
  19. _lapack_letter = ['s', 'd', 'c', 'z']
  20. _array_kind = {'i':0, 'l': 0, 'f': 0, 'd': 0, 'F': 1, 'D': 1}
  21. _array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
  22. _array_type = [['f', 'd'], ['F', 'D']]
  23.  
  24. def _commonType(*arrays):
  25.     kind = 0
  26. #    precision = 0
  27. #   force higher precision in lite version
  28.     precision = 1
  29.     for a in arrays:
  30.         t = a.typecode()
  31.         kind = max(kind, _array_kind[t])
  32.         precision = max(precision, _array_precision[t])
  33.     return _array_type[kind][precision]
  34.  
  35. def _castCopyAndTranspose(type, *arrays):
  36.     cast_arrays = ()
  37.     for a in arrays:
  38.         if a.typecode() == type:
  39.             cast_arrays = cast_arrays + (copy.copy(Numeric.transpose(a)),)
  40.         else:
  41.             cast_arrays = cast_arrays + (copy.copy(
  42.                                        Numeric.transpose(a).astype(type)),)
  43.     if len(cast_arrays) == 1:
  44.             return cast_arrays[0]
  45.     else:
  46.         return cast_arrays
  47.  
  48. # _fastCopyAndTranpose is an optimized version of _castCopyAndTranspose.
  49. # It assumes the input is 2D (as all the calls in here are).
  50.  
  51. _fastCT = multiarray._fastCopyAndTranspose
  52.  
  53. def _fastCopyAndTranspose(type, *arrays):
  54.     cast_arrays = ()
  55.     for a in arrays:
  56.         if a.typecode() == type:
  57.             cast_arrays = cast_arrays + (_fastCT(a),)
  58.         else:
  59.             cast_arrays = cast_arrays + (_fastCT(a.astype(type)),)
  60.     if len(cast_arrays) == 1:
  61.             return cast_arrays[0]
  62.     else:
  63.         return cast_arrays
  64.  
  65. def _assertRank2(*arrays):
  66.     for a in arrays:
  67.         if len(a.shape) != 2:
  68.             raise LinAlgError, 'Array must be two-dimensional'
  69.  
  70. def _assertSquareness(*arrays):
  71.     for a in arrays:
  72.         if max(a.shape) != min(a.shape):
  73.             raise LinAlgError, 'Array must be square'
  74.     
  75.  
  76. # Linear equations
  77.         
  78. def solve_linear_equations(a, b):
  79.     one_eq = len(b.shape) == 1
  80.     if one_eq:
  81.         b = b[:, Numeric.NewAxis]
  82.     _assertRank2(a, b)
  83.     _assertSquareness(a)
  84.     n_eq = a.shape[0]
  85.     n_rhs = b.shape[1]
  86.     if n_eq != b.shape[0]:
  87.         raise LinAlgError, 'Incompatible dimensions'
  88.     t =_commonType(a, b)
  89. #    lapack_routine = _findLapackRoutine('gesv', t)
  90.     if _array_kind[t] == 1: # Complex routines take different arguments
  91.         lapack_routine = lapack_lite.zgesv
  92.     else:
  93.         lapack_routine = lapack_lite.dgesv
  94.     a, b = _fastCopyAndTranspose(t, a, b)
  95.     pivots = Numeric.zeros(n_eq, 'i')
  96.     results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
  97.     if results['info'] > 0:
  98.         raise LinAlgError, 'Singular matrix'
  99.     if one_eq:
  100.         return Numeric.ravel(b) # I see no need to copy here
  101.     else:
  102.         return multiarray.transpose(b) # no need to copy
  103.  
  104.  
  105. # Matrix inversion
  106.  
  107. def inverse(a):
  108.     return solve_linear_equations(a, Numeric.identity(a.shape[0]))
  109.  
  110. # Cholesky decomposition
  111.  
  112. def cholesky_decomposition(a):
  113.     _assertRank2(a)
  114.     _assertSquareness(a)
  115.     t =_commonType(a)
  116.     a = _castCopyAndTranspose(t, a)
  117.     m = a.shape[0]
  118.     n = a.shape[1]
  119.     if _array_kind[t] == 1:
  120.         lapack_routine = lapack_lite.zpotrf
  121.     else:
  122.         lapack_routine = lapack_lite.dpotrf
  123.     results = lapack_routine('L', n, a, m, 0)
  124.     if results['info'] > 0:
  125.         raise LinAlgError, 'Matrix is not positive definite - Cholesky decomposition cannot be computed'
  126.     return copy.copy(Numeric.transpose(MLab.triu(a,k=0)))
  127.  
  128.  
  129. # Eigenvalues
  130.  
  131. def eigenvalues(a):
  132.     _assertRank2(a)
  133.     _assertSquareness(a)
  134.     t =_commonType(a)
  135.     real_t = _array_type[0][_array_precision[t]]
  136.     a = _fastCopyAndTranspose(t, a)
  137.     n = a.shape[0]
  138.     dummy = Numeric.zeros((1,), t)
  139.     if _array_kind[t] == 1: # Complex routines take different arguments
  140.         lapack_routine = lapack_lite.zgeev
  141.         w = Numeric.zeros((n,), t)
  142.         rwork = Numeric.zeros((n,),real_t)
  143.         lwork = 1
  144.         work = Numeric.zeros((lwork,), t)
  145.         results = lapack_routine('N', 'N', n, a, n, w,
  146.                                  dummy, 1, dummy, 1, work, -1, rwork, 0)
  147.         lwork = int(abs(work[0]))
  148.         work = Numeric.zeros((lwork,), t)
  149.         results = lapack_routine('N', 'N', n, a, n, w,
  150.                                  dummy, 1, dummy, 1, work, lwork, rwork, 0)
  151.     else:
  152.         lapack_routine = lapack_lite.dgeev
  153.         wr = Numeric.zeros((n,), t)
  154.         wi = Numeric.zeros((n,), t)
  155.         lwork = 1
  156.         work = Numeric.zeros((lwork,), t)
  157.         results = lapack_routine('N', 'N', n, a, n, wr, wi,
  158.                                  dummy, 1, dummy, 1, work, -1, 0)
  159.         lwork = int(work[0])
  160.         work = Numeric.zeros((lwork,), t)
  161.         results = lapack_routine('N', 'N', n, a, n, wr, wi,
  162.                                  dummy, 1, dummy, 1, work, lwork, 0)
  163.         if Numeric.logical_and.reduce(Numeric.equal(wi, 0.)):
  164.             w = wr
  165.         else:
  166.             w = wr+1j*wi
  167.     if results['info'] > 0:
  168.         raise LinAlgError, 'Eigenvalues did not converge'
  169.     return w
  170.  
  171.  
  172. def Heigenvalues(a, UPLO='L'):
  173.     _assertRank2(a)
  174.     _assertSquareness(a)
  175.     t =_commonType(a)
  176.     real_t = _array_type[0][_array_precision[t]]
  177.     a = _castCopyAndTranspose(t, a)
  178.     n = a.shape[0]
  179.     liwork = 5*n+3 
  180.     iwork = Numeric.zeros((liwork,),'i')
  181.     if _array_kind[t] == 1: # Complex routines take different arguments
  182.         lapack_routine = lapack_lite.zheevd
  183.         w = Numeric.zeros((n,), real_t)
  184.         lwork = 1
  185.         work = Numeric.zeros((lwork,), t)
  186.         lrwork = 1
  187.         rwork = Numeric.zeros((lrwork,),real_t)
  188.         results = lapack_routine('N', UPLO, n, a, n,w, work, -1, rwork, -1, iwork, liwork,  0)
  189.         lwork = int(abs(work[0]))
  190.         work = Numeric.zeros((lwork,), t)
  191.         lrwork = int(rwork[0])
  192.         rwork = Numeric.zeros((lrwork,),real_t)
  193.         results = lapack_routine('N', UPLO, n, a, n,w, work, lwork, rwork, lrwork, iwork, liwork,  0)
  194.     else:
  195.         lapack_routine = lapack_lite.dsyevd
  196.         w = Numeric.zeros((n,), t)
  197.         lwork = 1
  198.         work = Numeric.zeros((lwork,), t)
  199.         results = lapack_routine('N', UPLO, n, a, n,w, work, -1, iwork, liwork, 0)
  200.         lwork = int(work[0])
  201.         work = Numeric.zeros((lwork,), t)
  202.         results = lapack_routine('N', UPLO, n, a, n,w, work, lwork, iwork, liwork, 0)
  203.     if results['info'] > 0:
  204.         raise LinAlgError, 'Eigenvalues did not converge'
  205.     return w
  206.  
  207. # Eigenvectors
  208.  
  209. def eigenvectors(a):
  210.     """eigenvectors(a) returns u,v  where u is the eigenvalues and
  211. v is a matrix of eigenvectors with vector v[i] corresponds to 
  212. eigenvalue u[i].  Satisfies the equation dot(a, v[i]) = u[i]*v[i]
  213. """
  214.     _assertRank2(a)
  215.     _assertSquareness(a)
  216.     t =_commonType(a)
  217.     real_t = _array_type[0][_array_precision[t]]
  218.     a = _fastCopyAndTranspose(t, a)
  219.     n = a.shape[0]
  220.     dummy = Numeric.zeros((1,), t)
  221.     if _array_kind[t] == 1: # Complex routines take different arguments
  222.         lapack_routine = lapack_lite.zgeev
  223.         w = Numeric.zeros((n,), t)
  224.         v = Numeric.zeros((n,n), t)
  225.         lwork = 1
  226.         work = Numeric.zeros((lwork,),t)
  227.         rwork = Numeric.zeros((2*n,),real_t)
  228.         results = lapack_routine('N', 'V', n, a, n, w,
  229.                                   dummy, 1, v, n, work, -1, rwork, 0)
  230.         lwork = int(abs(work[0]))
  231.         work = Numeric.zeros((lwork,),t)
  232.         results = lapack_routine('N', 'V', n, a, n, w,
  233.                                   dummy, 1, v, n, work, lwork, rwork, 0)
  234.     else:
  235.         lapack_routine = lapack_lite.dgeev
  236.         wr = Numeric.zeros((n,), t)
  237.         wi = Numeric.zeros((n,), t)
  238.         vr = Numeric.zeros((n,n), t)
  239.         lwork = 1
  240.         work = Numeric.zeros((lwork,),t)
  241.         results = lapack_routine('N', 'V', n, a, n, wr, wi,
  242.                                   dummy, 1, vr, n, work, -1, 0)
  243.         lwork = int(work[0])
  244.         work = Numeric.zeros((lwork,),t)
  245.         results = lapack_routine('N', 'V', n, a, n, wr, wi,
  246.                                   dummy, 1, vr, n, work, lwork, 0)
  247.         if Numeric.logical_and.reduce(Numeric.equal(wi, 0.)):
  248.             w = wr
  249.             v = vr
  250.         else:
  251.             w = wr+1j*wi
  252.             v = Numeric.array(vr,Numeric.Complex)
  253.             ind = Numeric.nonzero(
  254.                           Numeric.equal(
  255.                               Numeric.equal(wi,0.0) # true for real e-vals
  256.                                        ,0)          # true for complex e-vals
  257.                                  )                  # indices of complex e-vals
  258.             for i in range(len(ind)/2):
  259.                 v[ind[2*i]] = vr[ind[2*i]] + 1j*vr[ind[2*i+1]]
  260.                 v[ind[2*i+1]] = vr[ind[2*i]] - 1j*vr[ind[2*i+1]]
  261.     if results['info'] > 0:
  262.         raise LinAlgError, 'Eigenvalues did not converge'
  263.     return w,v
  264.  
  265.  
  266. def Heigenvectors(a, UPLO='L'):
  267.     _assertRank2(a)
  268.     _assertSquareness(a)
  269.     t =_commonType(a)
  270.     real_t = _array_type[0][_array_precision[t]]
  271.     a = _castCopyAndTranspose(t, a)
  272.     n = a.shape[0]
  273.     liwork = 5*n+3
  274.     iwork = Numeric.zeros((liwork,),'i')
  275.     if _array_kind[t] == 1: # Complex routines take different arguments
  276.         lapack_routine = lapack_lite.zheevd
  277.         w = Numeric.zeros((n,), real_t)
  278.         lwork = 1
  279.         work = Numeric.zeros((lwork,), t)
  280.         lrwork = 1
  281.         rwork = Numeric.zeros((lrwork,),real_t)
  282.         results = lapack_routine('V', UPLO, n, a, n,w, work, -1, rwork, -1, iwork, liwork,  0)
  283.         lwork = int(abs(work[0]))
  284.         work = Numeric.zeros((lwork,), t)
  285.         lrwork = int(rwork[0])
  286.         rwork = Numeric.zeros((lrwork,),real_t)
  287.         results = lapack_routine('V', UPLO, n, a, n,w, work, lwork, rwork, lrwork, iwork, liwork,  0)
  288.     else:
  289.         lapack_routine = lapack_lite.dsyevd
  290.         w = Numeric.zeros((n,), t)
  291.         lwork = 1
  292.         work = Numeric.zeros((lwork,),t)
  293.         results = lapack_routine('V', UPLO, n, a, n,w, work, -1, iwork, liwork, 0)
  294.         lwork = int(work[0])
  295.         work = Numeric.zeros((lwork,),t)
  296.         results = lapack_routine('V', UPLO, n, a, n,w, work, lwork, iwork, liwork, 0)
  297.     if results['info'] > 0:
  298.         raise LinAlgError, 'Eigenvalues did not converge'
  299.     return (w,a)
  300.  
  301.    
  302. # Singular value decomposition
  303.  
  304. def singular_value_decomposition(a, full_matrices = 0):
  305.     _assertRank2(a)
  306.     n = a.shape[1]
  307.     m = a.shape[0]
  308.     t =_commonType(a)
  309.     real_t = _array_type[0][_array_precision[t]]
  310.     a = _fastCopyAndTranspose(t, a)
  311.     if full_matrices:
  312.         nu = m
  313.         nvt = n
  314.         option = 'A'
  315.     else:
  316.         nu = min(n,m)
  317.         nvt = min(n,m)
  318.         option = 'S'
  319.     s = Numeric.zeros((min(n,m),), real_t)
  320.     u = Numeric.zeros((nu, m), t)
  321.     vt = Numeric.zeros((n, nvt), t)
  322.     iwork = Numeric.zeros((8*min(m,n),), 'i')
  323.     if _array_kind[t] == 1: # Complex routines take different arguments
  324.         lapack_routine = lapack_lite.zgesdd
  325.         rwork = Numeric.zeros((5*min(m,n)*min(m,n) + 5*min(m,n),), real_t)
  326.         lwork = 1
  327.         work = Numeric.zeros((lwork,), t)
  328.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt,
  329.                                  work, -1, rwork, iwork, 0)
  330.         lwork = int(abs(work[0]))
  331.         work = Numeric.zeros((lwork,), t)
  332.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt,
  333.                                  work, lwork, rwork, iwork, 0)
  334.     else:
  335.         lapack_routine = lapack_lite.dgesdd
  336.         lwork = 1
  337.         work = Numeric.zeros((lwork,), t)
  338.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt,
  339.                                  work, -1, iwork, 0)
  340.         lwork = int(work[0])
  341.         work = Numeric.zeros((lwork,), t)
  342.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt,
  343.                                  work, lwork, iwork, 0)
  344.     if results['info'] > 0:
  345.         raise LinAlgError, 'SVD did not converge'
  346.     return multiarray.transpose(u), s, multiarray.transpose(vt) # why copy here?
  347.  
  348.  
  349. # Generalized inverse
  350.  
  351. def generalized_inverse(a, rcond = 1.e-10):
  352.     a = Numeric.array(a, copy=0)
  353.     if a.typecode() in Numeric.typecodes['Complex']:
  354.         a = Numeric.conjugate(a)
  355.     u, s, vt = singular_value_decomposition(a, 0)
  356.     m = u.shape[0]
  357.     n = vt.shape[1]
  358.     cutoff = rcond*Numeric.maximum.reduce(s)
  359.     for i in range(min(n,m)):
  360.         if s[i] > cutoff:
  361.             s[i] = 1./s[i]
  362.         else:
  363.             s[i] = 0.;
  364.     return Numeric.dot(Numeric.transpose(vt),
  365.                        s[:, Numeric.NewAxis]*Numeric.transpose(u))
  366.  
  367. # Determinant
  368.  
  369. def determinant(a):
  370.     _assertRank2(a)
  371.     _assertSquareness(a)
  372.     t =_commonType(a)
  373.     a = _fastCopyAndTranspose(t, a)
  374.     n = a.shape[0]
  375.     if _array_kind[t] == 1:
  376.         lapack_routine = lapack_lite.zgetrf
  377.     else:
  378.         lapack_routine = lapack_lite.dgetrf
  379.     pivots = Numeric.zeros((n,), 'i')
  380.     results = lapack_routine(n, n, a, n, pivots, 0)
  381.     sign = Numeric.add.reduce(Numeric.not_equal(pivots,
  382.                                                 Numeric.arrayrange(1, n+1))) % 2
  383.     return (1.-2.*sign)*Numeric.multiply.reduce(Numeric.diagonal(a))
  384.  
  385. # Linear Least Squares 
  386.         
  387. def linear_least_squares(a, b, rcond=1.e-10):
  388.     """solveLinearLeastSquares(a,b) returns x,resids,rank,s 
  389. where x minimizes 2-norm(|b - Ax|) 
  390.       resids is the sum square residuals
  391.       rank is the rank of A
  392.       s is an rank of the singual values of A in desending order
  393.  
  394. If b is a matrix then x is also a matrix with corresponding columns.
  395. If the rank of A is less than the number of columns of A or greater than
  396. the numer of rows, then residuals will be returned as an empty array
  397. otherwise resids = sum((b-dot(A,x)**2).
  398. Singular values less than s[0]*rcond are treated as zero.
  399. """
  400.     one_eq = len(b.shape) == 1
  401.     if one_eq:
  402.         b = b[:, Numeric.NewAxis]
  403.     _assertRank2(a, b)
  404.     m  = a.shape[0]
  405.     n  = a.shape[1]
  406.     n_rhs = b.shape[1]
  407.     ldb = max(n,m)
  408.     if m != b.shape[0]:
  409.         raise LinAlgError, 'Incompatible dimensions'
  410.     t =_commonType(a, b)
  411.     real_t = _array_type[0][_array_precision[t]]
  412.     bstar = Numeric.zeros((ldb,n_rhs),t)
  413.     bstar[:b.shape[0],:n_rhs] = copy.copy(b)
  414.     a,bstar = _castCopyAndTranspose(t, a, bstar)
  415.     s = Numeric.zeros((min(m,n),),real_t)
  416.     nlvl = max( 0, int( math.log( float(min( m,n ))/2. ) ) + 1 )
  417.     iwork = Numeric.zeros((3*min(m,n)*nlvl+11*min(m,n),), 'i')
  418.     if _array_kind[t] == 1: # Complex routines take different arguments
  419.         lapack_routine = lapack_lite.zgelsd
  420.         lwork = 1
  421.         rwork = Numeric.zeros((lwork,), real_t)
  422.         work = Numeric.zeros((lwork,),t)
  423.         results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
  424.                         0,work,-1,rwork,iwork,0 )
  425.         lwork = int(abs(work[0]))
  426.         rwork = Numeric.zeros((lwork,),real_t)
  427.         a_real = Numeric.zeros((m,n),real_t)
  428.         bstar_real = Numeric.zeros((ldb,n_rhs,),real_t)
  429.         results = lapack_lite.dgelsd( m, n, n_rhs, a_real, m, bstar_real,ldb , s, rcond,
  430.                         0,rwork,-1,iwork,0 )
  431.         lrwork = int(rwork[0])
  432.         work = Numeric.zeros((lwork,), t)
  433.         rwork = Numeric.zeros((lrwork,), real_t)
  434.         results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
  435.                         0,work,lwork,rwork,iwork,0 )
  436.     else:
  437.         lapack_routine = lapack_lite.dgelsd
  438.         lwork = 1
  439.         work = Numeric.zeros((lwork,), t)
  440.         results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
  441.                         0,work,-1,iwork,0 )
  442.         lwork = int(work[0])
  443.         work = Numeric.zeros((lwork,), t)
  444.         results = lapack_routine( m, n, n_rhs, a, m, bstar,ldb , s, rcond,
  445.                         0,work,lwork,iwork,0 )
  446.     if results['info'] > 0:
  447.         raise LinAlgError, 'SVD did not converge in Linear Least Squares'
  448.     resids = Numeric.array([],t)
  449.     if one_eq:
  450.         x = copy.copy(Numeric.ravel(bstar)[:n])
  451.         if (results['rank']==n) and (m>n):
  452.             resids = Numeric.array([Numeric.sum((Numeric.ravel(bstar)[n:])**2)])
  453.     else:
  454.         x = copy.copy(Numeric.transpose(bstar)[:n,:])
  455.         if (results['rank']==n) and (m>n):
  456.             resids = copy.copy(Numeric.sum((Numeric.transpose(bstar)[n:,:])**2))
  457.     return x,resids,results['rank'],copy.copy(s[:min(n,m)]) 
  458.  
  459.  
  460. if __name__ == '__main__':
  461.         from Numeric import *
  462.  
  463.         def test(a, b):
  464.  
  465.                 print "All numbers printed should be (almost) zero:"
  466.  
  467.                 x = solve_linear_equations(a, b)
  468.                 check = b - matrixmultiply(a, x)
  469.                 print check
  470.  
  471.  
  472.                 a_inv = inverse(a)
  473.                 check = matrixmultiply(a, a_inv)-identity(a.shape[0])
  474.                 print check
  475.  
  476.  
  477.                 ev = eigenvalues(a)
  478.  
  479.                 evalues, evectors = eigenvectors(a)
  480.                 check = ev-evalues
  481.                 print check
  482.  
  483.                 evectors = transpose(evectors)
  484.                 check = matrixmultiply(a, evectors)-evectors*evalues
  485.                 print check
  486.  
  487.  
  488.                 u, s, vt = singular_value_decomposition(a)
  489.                 check = a - Numeric.matrixmultiply(u*s, vt)
  490.                 print check
  491.  
  492.  
  493.                 a_ginv = generalized_inverse(a)
  494.                 check = matrixmultiply(a, a_ginv)-identity(a.shape[0])
  495.                 print check
  496.  
  497.  
  498.                 det = determinant(a)
  499.                 check = det-multiply.reduce(evalues)
  500.                 print check
  501.  
  502.                 x, residuals, rank, sv = linear_least_squares(a, b)
  503.                 check = b - matrixmultiply(a, x)
  504.                 print check
  505.                 print rank-a.shape[0]
  506.                 print sv-s
  507.  
  508.         a = array([[1.,2.], [3.,4.]])
  509.         b = array([2., 1.])
  510.         test(a, b)
  511.  
  512.         a = a+0j
  513.         b = b+0j
  514.         test(a, b)
  515.  
  516.